Eye-body allometry in each species

We first examine static post-metamorphosis eye diameter-body length allometry across 16 species of deep-sea sergestids. Allometric relationships were fit using standardised major axis regression via the smatr v.3.4.8 package. For each species, standard model checks, model fits, and an interactive plot of data and fits are shown. Hover over the plots to see individual data values. Points are colored by the source of the sample (fresh collection or preserved in ethanol or paraformaldehyde).

Note that for 3 of the species, a large/unreasonable outlier was present. We identified these by eye and they have been removed from the data subsets and analyses, but the outlier value is noted under each species header for reference, and a second plot including the outlier as a red diamond is included.

#import specimen data
specimens <- data.frame(read.csv("../Data/cleaned data/specimen_data_tidy.csv", header=TRUE))
#set colors for fresh and preserved samples
col_pres <- c("ethanol" = "#a6cee3",
              "fresh" = "#b2df8a",
              "paraformaldehyde" = "#1f78b4")

Allosergestes pectinatus

# Subset data -----
apec <- specimens %>% filter(genus_species == "Allosergestes_pectinatus")

# Fit SMA model ------
sma_apec <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = apec, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_apec, which = "default",type = "o") 
plot(sma_apec, which = "residual",type = "o")
plot(sma_apec, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_apec)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = apec, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
## 
## H0 : variables uncorrelated
## R-squared : 0.8156794 
## P-value : 9.643e-06
#save coefficients of fits as object
cc_sma_apec <- data.frame(coef(sma_apec))

# Make plot ------
plot_apec <- ggplot(apec, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Allosergestes pectinatus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_apec, aes(intercept = cc_sma_apec[1,1], slope = cc_sma_apec[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_apec[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_apec[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_apec)

Allosergestes sargassi

Note: this species has a big outlier that was removed for analyses and plots (eye 0.7, body 13.0).

# Subset data -----
#note: outlier removed
asar <- specimens %>% 
  filter(genus_species == "Allosergestes_sargassi") %>%
  filter(!(Eye_Diameter == round(0.67, 1) & Body_Length == round(12.97,1)))

# Fit SMA model ------
sma_asar <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = asar, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_asar, which = "default",type = "o") 
plot(sma_asar, which = "residual",type = "o")
plot(sma_asar, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_asar)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = asar, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
## 
## H0 : variables uncorrelated
## R-squared : 0.4037707 
## P-value : 2.3852e-05
#save coefficients of fits as object
cc_sma_asar <- data.frame(coef(sma_asar))

# Make plot ------
plot_asar <- ggplot(asar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Allosergestes sargassi") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_asar, aes(intercept = cc_sma_asar[1,1], slope = cc_sma_asar[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_asar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_asar[1,1], digits = 2), sep = " "), size = 3, col = "black")

#plot with outlier shown
plot_asar_out <- plot_asar +
  geom_point(aes(y=round(0.67, 1), x=round(12.97,1)), colour="red", pch = 18) +
  ggtitle("Allosergestes sargassi + outlier")
  
#plot
ggplotly(plot_asar)
#plot with outlier
ggplotly(plot_asar_out)

Challengerosergia hansjacobi

# Subset data -----
chan <- specimens %>% filter(genus_species == "Challengerosergia_hansjacobi")

# Fit SMA model ------
sma_chan <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = chan, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_chan, which = "default",type = "o") 
plot(sma_chan, which = "residual",type = "o")
plot(sma_chan, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_chan)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = chan, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
## 
## H0 : variables uncorrelated
## R-squared : 0.4129122 
## P-value : 3.3885e-06
#save coefficients of fits as object
cc_sma_chan <- data.frame(coef(sma_chan))

# Make plot ------
plot_chan <- ggplot(chan, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Challengerosergia hansjacobi") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_chan, aes(intercept = cc_sma_chan[1,1], slope = cc_sma_chan[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_chan[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_chan[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_chan)

Challengerosergia talismani

# Subset data -----
ctal <- specimens %>% filter(genus_species == "Challengerosergia_talismani")

# Fit SMA model ------
sma_ctal <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = ctal, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_ctal, which = "default",type = "o") 
plot(sma_ctal, which = "residual",type = "o")
plot(sma_ctal, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_ctal)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = ctal, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
## 
## H0 : variables uncorrelated
## R-squared : 0.8120894 
## P-value : 2.0471e-12
#save coefficients of fits as object
cc_sma_ctal <- data.frame(coef(sma_ctal))

# Make plot ------
plot_ctal <- ggplot(ctal, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Challengerosergia talismani") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_ctal, aes(intercept = cc_sma_ctal[1,1], slope = cc_sma_ctal[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_ctal[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_ctal[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_ctal)

Deosergestes corniculum

# Subset data -----
dcor <- specimens %>% filter(genus_species == "Deosergestes_corniculum")

# Fit SMA model ------
sma_dcor <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = dcor, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dcor, which = "default",type = "o") 
plot(sma_dcor, which = "residual",type = "o")
plot(sma_dcor, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_dcor)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dcor, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
## 
## H0 : variables uncorrelated
## R-squared : 0.8421211 
## P-value : 8.1768e-08
#save coefficients of fits as object
cc_sma_dcor <- data.frame(coef(sma_dcor))

# Make plot ------
plot_dcor <- ggplot(dcor, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Deosergestes corniculum") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_dcor, aes(intercept = cc_sma_dcor[1,1], slope = cc_sma_dcor[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dcor[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dcor[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_dcor)

Deosergestes henseni

# Subset data -----
dhen <- specimens %>% filter(genus_species == "Deosergestes_henseni")

# Fit SMA model ------
sma_dhen <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = dhen, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dhen, which = "default",type = "o") 
plot(sma_dhen, which = "residual",type = "o")
plot(sma_dhen, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_dhen)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dhen, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
## 
## H0 : variables uncorrelated
## R-squared : 0.3380975 
## P-value : 1.6982e-06
#save coefficients of fits as object
cc_sma_dhen <- data.frame(coef(sma_dhen))

# Make plot ------
plot_dhen <- ggplot(dhen, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Deosergestes henseni") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_dhen, aes(intercept = cc_sma_dhen[1,1], slope = cc_sma_dhen[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dhen[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dhen[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_dhen)

Eusergestes arcticus

# Subset data -----
euar <- specimens %>% filter(genus_species == "Eusergestes_arcticus")

# Fit SMA model ------
sma_euar <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = euar, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_euar, which = "default",type = "o") 
plot(sma_euar, which = "residual",type = "o")
plot(sma_euar, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_euar)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = euar, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
## 
## H0 : variables uncorrelated
## R-squared : 0.8368056 
## P-value : 0.029484
#save coefficients of fits as object
cc_sma_euar <- data.frame(coef(sma_euar))

# Make plot ------
plot_euar <- ggplot(euar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Eusergestes arcticus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_euar, aes(intercept = cc_sma_euar[1,1], slope = cc_sma_euar[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_euar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_euar[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_euar)

Gardinerosergia splendens

Note: this species has a big outlier that was removed for analyses and plots (eye 1.7, body 22.0).

# Subset data -----
gspl <- specimens %>% 
  filter(genus_species == "Gardinerosergia_splendens") %>%
  filter(!(Eye_Diameter == round(1.71, 1) & Body_Length == round(21.98,1)))

# Fit SMA model ------
sma_gspl <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = gspl, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_gspl, which = "default",type = "o") 
plot(sma_gspl, which = "residual",type = "o")
plot(sma_gspl, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_gspl)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = gspl, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
## 
## H0 : variables uncorrelated
## R-squared : 0.4597308 
## P-value : 1.1132e-06
#save coefficients of fits as object
cc_sma_gspl <- data.frame(coef(sma_gspl))

# Make plot ------
plot_gspl <- ggplot(gspl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Gardinerosergia splendens") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_gspl, aes(intercept = cc_sma_gspl[1,1], slope = cc_sma_gspl[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_gspl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_gspl[1,1], digits = 2), sep = " "), size = 3, col = "black")

#plot with outlier shown
plot_gspl_out <- plot_gspl +
  geom_point(aes(y=round(1.71,1), x=round(21.98,1)), colour="red", pch = 18) +
  ggtitle("Gardinerosergia splendens + outlier")
  
#plot
ggplotly(plot_gspl)
#plot with outlier
ggplotly(plot_gspl_out)

Neosergestes edwardsii

Note: this species has a big outlier that was removed for analyses and plots (eye 0.9, body 0.7). Also note that for this species, sample size was insufficient to estimate eye-body allometry (N = 6, p = 0.57 for correlation between eye size and body size).

# Subset data -----
nedw <- specimens %>% 
  filter(genus_species == "Neosergestes_edwardsii") %>%
  filter(!(Eye_Diameter == round(0.85,1) & Body_Length == round(0.72,1)))

# Fit SMA model ------
sma_nedw <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = nedw, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_nedw, which = "default",type = "o") 
plot(sma_nedw, which = "residual",type = "o")
plot(sma_nedw, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_nedw)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = nedw, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation      slope
## estimate     1.4244363 -1.2712830
## lower limit -0.7533312 -3.8007834
## upper limit  3.6022037 -0.4252177
## 
## H0 : variables uncorrelated
## R-squared : 0.08540156 
## P-value : 0.57413
#save coefficients of fits as object
cc_sma_nedw <- data.frame(coef(sma_nedw))

# Make plot ------
plot_nedw <- ggplot(nedw, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Neosergestes edwardsii") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) 

#plot with outlier shown
plot_nedw_out <- plot_nedw +
  geom_point(aes(x=round(0.85,1), y=round(0.72,1)), colour="red", pch = 18) +
  ggtitle("Neosergestes edwardsii + outlier")
  
#plot
ggplotly(plot_nedw)
#plot with outlier
ggplotly(plot_nedw_out)

Parasergestes armatus

# Subset data -----
parm <- specimens %>% filter(genus_species == "Parasergestes_armatus")

# Fit SMA model ------
sma_parm <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = parm, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_parm, which = "default",type = "o") 
plot(sma_parm, which = "residual",type = "o")
plot(sma_parm, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_parm)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = parm, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
## 
## H0 : variables uncorrelated
## R-squared : 0.5670657 
## P-value : 1.5744e-06
#save coefficients of fits as object
cc_sma_parm <- data.frame(coef(sma_parm))

# Make plot ------
plot_parm <- ggplot(parm, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Parasergestes armatus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_parm, aes(intercept = cc_sma_parm[1,1], slope = cc_sma_parm[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_parm[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_parm[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_parm)

Parasergestes vigilax

# Subset data -----
pvig <- specimens %>% filter(genus_species == "Parasergestes_vigilax")

# Fit SMA model ------
sma_pvig <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = pvig, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pvig, which = "default",type = "o") 
plot(sma_pvig, which = "residual",type = "o")
plot(sma_pvig, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_pvig)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pvig, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
## 
## H0 : variables uncorrelated
## R-squared : 0.5274941 
## P-value : 0.00096145
#save coefficients of fits as object
cc_sma_pvig <- data.frame(coef(sma_pvig))

# Make plot ------
plot_pvig <- ggplot(pvig, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Parasergestes vigilax") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_pvig, aes(intercept = cc_sma_pvig[1,1], slope = cc_sma_pvig[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pvig[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pvig[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_pvig)

Phorcosergia grandis

# Subset data -----
pgra <- specimens %>% filter(genus_species == "Phorcosergia_grandis")

# Fit SMA model ------
sma_pgra <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = pgra, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pgra, which = "default",type = "o") 
plot(sma_pgra, which = "residual",type = "o")
plot(sma_pgra, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_pgra)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pgra, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
## 
## H0 : variables uncorrelated
## R-squared : 0.8798127 
## P-value : < 2.22e-16
#save coefficients of fits as object
cc_sma_pgra <- data.frame(coef(sma_pgra))

# Make plot ------
plot_pgra <- ggplot(pgra, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Phorcosergia grandis") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_pgra, aes(intercept = cc_sma_pgra[1,1], slope = cc_sma_pgra[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pgra[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pgra[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_pgra)

Robustosergia robusta

# Subset data -----
rrob <- specimens %>% filter(genus_species == "Robustosergia_robusta")

# Fit SMA model ------
sma_rrob <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = rrob, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rrob, which = "default",type = "o") 
plot(sma_rrob, which = "residual",type = "o")
plot(sma_rrob, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_rrob)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rrob, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
## 
## H0 : variables uncorrelated
## R-squared : 0.6989428 
## P-value : 6.7833e-07
#save coefficients of fits as object
cc_sma_rrob <- data.frame(coef(sma_rrob))

# Make plot ------
plot_rrob <- ggplot(rrob, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Robustosergia robusta") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_rrob, aes(intercept = cc_sma_rrob[1,1], slope = cc_sma_rrob[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rrob[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rrob[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_rrob)

Robustosergia regalis

# Subset data -----
rreg <- specimens %>% filter(genus_species == "Robustosergia_regalis")

# Fit SMA model ------
sma_rreg <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = rreg, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rreg, which = "default",type = "o") 
plot(sma_rreg, which = "residual",type = "o")
plot(sma_rreg, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_rreg)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rreg, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
## 
## H0 : variables uncorrelated
## R-squared : 0.8280327 
## P-value : 2.4971e-15
#save coefficients of fits as object
cc_sma_rreg <- data.frame(coef(sma_rreg))

# Make plot ------
plot_rreg <- ggplot(rreg, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Robustosergia regalis") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_rreg, aes(intercept = cc_sma_rreg[1,1], slope = cc_sma_rreg[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rreg[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rreg[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_rreg)

Sergestes atlanticus

# Subset data -----
satl <- specimens %>% filter(genus_species == "Sergestes_atlanticus")

# Fit SMA model ------
sma_satl <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = satl, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_satl, which = "default",type = "o") 
plot(sma_satl, which = "residual",type = "o")
plot(sma_satl, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_satl)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = satl, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
## 
## H0 : variables uncorrelated
## R-squared : 0.7067609 
## P-value : 8.6609e-05
#save coefficients of fits as object
cc_sma_satl <- data.frame(coef(sma_satl))

# Make plot ------
plot_satl <- ggplot(satl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Sergestes atlanticus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_satl, aes(intercept = cc_sma_satl[1,1], slope = cc_sma_satl[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_satl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_satl[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_satl)

Sergia tenuiremis

# Subset data -----
sten <- specimens %>% filter(genus_species == "Sergia_tenuiremis")

# Fit SMA model ------
sma_sten <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = sten, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_sten, which = "default",type = "o") 
plot(sma_sten, which = "residual",type = "o")
plot(sma_sten, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_sten)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = sten, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
## 
## H0 : variables uncorrelated
## R-squared : 0.2880646 
## P-value : 0.0032339
#save coefficients of fits as object
cc_sma_sten <- data.frame(coef(sma_sten))

# Make plot ------
plot_sten <- ggplot(sten, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Sergia tenuiremis") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_sten, aes(intercept = cc_sma_sten[1,1], slope = cc_sma_sten[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_sten[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_sten[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_sten)

Figure for static eye-body scaling across sergestid species

#make figure panels
fig.a <- plot_dcor + theme(legend.position = "none", axis.title.x = element_blank())
fig.b <- plot_dhen + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.c <- plot_apec + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.d <- plot_asar + theme(legend.position = "none", axis.title.x = element_blank())
fig.e <- plot_satl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.f <- plot_pvig + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.g <- plot_parm + theme(legend.position = "none", axis.title.x = element_blank())
fig.h <- plot_euar + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.i <- plot_gspl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.j <- plot_rreg + theme(legend.position = "none", axis.title.x = element_blank())
fig.k <- plot_rrob + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.l <- plot_pgra + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.m <- plot_sten + theme(legend.position = "none")
fig.n <- plot_ctal  + theme(legend.position = "none", axis.title.y = element_blank())
fig.o <- plot_chan + theme(legend.position = "none", axis.title.y = element_blank())

# arrange panels in figure
plots <- plot_grid(fig.a, fig.b, fig.c, fig.d, fig.e, fig.f, fig.g, fig.h, fig.i, fig.j, fig.k, fig.l, fig.m, fig.n, fig.o, #list of plots to arrange in grid
           align = 'vh', #align horizontally and vertically
           axis = 'lb',
           labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I","J", "K", "L", "M", "N", "O"), #panel labels for figure
           hjust = -1, #adjustment for panel labels
           nrow = 5) #number of rows in grids

#extract legend
leg <- get_legend(fig.e + theme(legend.position="bottom"))


#view figure
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .01), align = 'vh')

#export figure
#png("../Figures/dev-allometry.png", width = 10, height = 12, units = "in", res = 500)
pdf("../Figures/Fig-S1.pdf", width = 10, height = 12)
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .05), align = 'vh')
dev.off()

Test whether species differ in static allometry

We’ve looked at how eyes scale with body size within each species, but are there significant differences in static allometry across species? Here we use two approaches: first, we compare linear mixed models to see which best describes our data. Then, we look at pairwise comparisons of SMA regressions among species to see which are driving these differences.

Linear mixed model to test whether slopes vary across species

Here, we use linear mixed models in lme4 v.1.1.25 to test whether species show differences in static eye-body allometric slopes. We fit 1) a constant-slope, variable intercept model to our eye and body size data with species as a random effect and 2) a variable slope model with the same effects, and then compared model fits by AIC scores. Both models are fit using reduced maximum likelihood (REML), as this approach is better than maximum likelihood (ML) for detecting differences in random effects among models, which is what we are testing here (ML is better for fixed effects).

First we removed the three outliers and the species with insufficient sampling (N. edwardsii) from our full specimen dataset, then ran each model.

#Remove outliers from specimen dataset
specimens_ed <- specimens %>%
  filter(!(genus_species == "Allosergestes_sargassi" & Eye_Diameter == 0.7 & Body_Length == 13.0)) %>% #(0.67, 12.97 raw)
  filter(!(genus_species == "Gardinerosergia_splendens" & Eye_Diameter == 1.7 & Body_Length == 22.0)) %>% #(1.71, 21.98 raw)
  filter(!(genus_species == "Neosergestes_edwardsii" & Eye_Diameter == 0.8 & Body_Length == 0.7)) #(0.85, 0.72 raw)

#Export data without outliers
write.csv(specimens_ed, file = "../Data/cleaned data/specimen_data_tidy_no-outliers.csv", row.names = FALSE)

#Remove species with insufficient sampling
specimens_test <- specimens_ed %>%
  filter(genus_species != "Neosergestes_edwardsii") %>%
  mutate(species_text = as.factor(paste(Genus, Species, sep = " "))) # Add labels for species text

Fixed slope (variable intercepts) model

# variable intercepts model (species can have different intercepts but share mean slope) -------
lme_fixed <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1|genus_species), 
                  data = specimens_test)

# summary 
summary(lme_fixed)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
##    Data: specimens_test
## 
## REML criterion at convergence: -1091.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6547 -0.5660  0.0425  0.5699  3.7587 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  genus_species (Intercept) 0.004257 0.06525 
##  Residual                  0.004549 0.06745 
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        -1.24683    0.05121  -24.35
## log10(Body_Length)  0.80642    0.03122   25.83
## 
## Correlation of Fixed Effects:
##             (Intr)
## lg10(Bdy_L) -0.941
# intercepts (variable) and slope (constant)
coef(lme_fixed)
## $genus_species
##                              (Intercept) log10(Body_Length)
## Allosergestes_pectinatus       -1.337988          0.8064169
## Allosergestes_sargassi         -1.295962          0.8064169
## Challengerosergia_hansjacobi   -1.164606          0.8064169
## Challengerosergia_talismani    -1.216875          0.8064169
## Deosergestes_corniculum        -1.328655          0.8064169
## Deosergestes_henseni           -1.278582          0.8064169
## Eusergestes_arcticus           -1.246866          0.8064169
## Gardinerosergia_splendens      -1.158266          0.8064169
## Parasergestes_armatus          -1.325352          0.8064169
## Parasergestes_vigilax          -1.299386          0.8064169
## Phorcosergia_grandis           -1.246064          0.8064169
## Robustosergia_regalis          -1.176361          0.8064169
## Robustosergia_robusta          -1.155267          0.8064169
## Sergestes_atlanticus           -1.234190          0.8064169
## Sergia_tenuiremis              -1.238014          0.8064169
## 
## attr(,"class")
## [1] "coef.mer"
# AIC
AIC(lme_fixed)
## [1] -1083.888
#fixed effects
fixef(lme_fixed)
##        (Intercept) log10(Body_Length) 
##         -1.2468288          0.8064169
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_fixed, level = 0.95)
##                          2.5 %      97.5 %
## .sig01              0.04428477  0.09578174
## .sigma              0.06313279  0.07211335
## (Intercept)        -1.34880013 -1.14643360
## log10(Body_Length)  0.74536494  0.86929119

Variable slopes model

# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_variable <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species), 
                     data = specimens_test)

# summary
summary(lme_variable)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |  
##     genus_species)
##    Data: specimens_test
## 
## REML criterion at convergence: -1125.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4279 -0.5359  0.0306  0.5910  4.2687 
## 
## Random effects:
##  Groups        Name               Variance Std.Dev. Corr 
##  genus_species (Intercept)        0.232807 0.48250       
##                log10(Body_Length) 0.104532 0.32331  -0.99
##  Residual                         0.003977 0.06306       
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        -1.35693    0.13739  -9.876
## log10(Body_Length)  0.89778    0.09181   9.778
## 
## Correlation of Fixed Effects:
##             (Intr)
## lg10(Bdy_L) -0.993
# intercepts (variable) and slope (constant)
coef(lme_variable)
## $genus_species
##                              (Intercept) log10(Body_Length)
## Allosergestes_pectinatus      -2.2120964          1.4762026
## Allosergestes_sargassi        -1.7513868          1.1366814
## Challengerosergia_hansjacobi  -1.0760213          0.7484245
## Challengerosergia_talismani   -1.4111004          0.9318641
## Deosergestes_corniculum       -1.1953553          0.7273404
## Deosergestes_henseni          -0.7843247          0.4963146
## Eusergestes_arcticus          -1.2644263          0.8196651
## Gardinerosergia_splendens     -0.8716955          0.6169987
## Parasergestes_armatus         -1.3258605          0.8075505
## Parasergestes_vigilax         -2.0269469          1.3719045
## Phorcosergia_grandis          -1.3536721          0.8693323
## Robustosergia_regalis         -1.2328227          0.8412590
## Robustosergia_robusta         -1.1498050          0.8037740
## Sergestes_atlanticus          -1.8945693          1.2766884
## Sergia_tenuiremis             -0.8038697          0.5427274
## 
## attr(,"class")
## [1] "coef.mer"
# AIC
AIC(lme_variable)
## [1] -1113.879
#fixed effects
fixef(lme_variable)
##        (Intercept) log10(Body_Length) 
##         -1.3569302          0.8977818
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_variable, level = 0.95)
##                          2.5 %      97.5 %
## .sig01              0.27582256  0.74974111
## .sig02             -0.99580186 -0.97348803
## .sig03              0.17701121  0.50987554
## .sigma              0.05903897  0.06767603
## (Intercept)        -1.64344107 -1.07983292
## log10(Body_Length)  0.71455352  1.09409513

Note that the parameter estimates for slopes and intercepts for species groups are a bit different than those we estimated via OLS. That’s because here, the model is looking across all the data for all the species, and finding a common mean slope/intercept as well as the group effects.

Compare models

Then we can compare the two models via AIC and log likelihood.

#model AIC comparison
AIC(lme_fixed, lme_variable)
##              df       AIC
## lme_fixed     4 -1083.888
## lme_variable  6 -1113.879
#log likelihood comparison
anova(lme_fixed, lme_variable, refit=FALSE)
## Data: specimens_test
## Models:
## lme_fixed: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
## lme_variable: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species)
##              npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## lme_fixed       4 -1083.9 -1067.5 545.94  -1091.9                         
## lme_variable    6 -1113.9 -1089.2 562.94  -1125.9 33.991  2  4.159e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here we find by comparing AIC scores and log likelihoods that the model with variable slopes is better supported than the one with fixed slopes. We have evidence here that species differ significantly in the slopes of their static eye-body allometries rather than following a shared growth trajectory.

Test for significant effect of preservation as a covariate

Here, we add preservation method as a fixed effect into our linear mixed models of eye size vs. body size with species as a random effect. We have already selected the variable slopes model, so here we use a maximum likelihood approach (ML) to fit models that 1) don’t include preservation as a covariate and 2) do include preservation as a covariate and then compare model fits. We include preservation as a fixed rather than a random effect because random effects should have more than 5 levels, and preservation type has only 3.

Model without preservation type

# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_nopres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species), 
                   data = specimens_test,
                   REML = FALSE) #uses ML instead of REML

# summary
summary(lme_nopres)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |  
##     genus_species)
##    Data: specimens_test
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1123.3  -1098.6    567.6  -1135.3      444 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4474 -0.5397  0.0290  0.5871  4.2627 
## 
## Random effects:
##  Groups        Name               Variance Std.Dev. Corr 
##  genus_species (Intercept)        0.205885 0.45375       
##                log10(Body_Length) 0.091871 0.30310  -0.99
##  Residual                         0.003983 0.06311       
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        -1.35193    0.13036  -10.37
## log10(Body_Length)  0.89371    0.08686   10.29
## 
## Correlation of Fixed Effects:
##             (Intr)
## lg10(Bdy_L) -0.993
# intercepts (variable) and slope (constant)
coef(lme_nopres)
## $genus_species
##                              (Intercept) log10(Body_Length)
## Allosergestes_pectinatus      -2.1804353          1.4518078
## Allosergestes_sargassi        -1.7432761          1.1308509
## Challengerosergia_hansjacobi  -1.0817350          0.7520768
## Challengerosergia_talismani   -1.4101465          0.9312293
## Deosergestes_corniculum       -1.1915651          0.7253703
## Deosergestes_henseni          -0.7877696          0.4985142
## Eusergestes_arcticus          -1.2643317          0.8197744
## Gardinerosergia_splendens     -0.8765164          0.6200697
## Parasergestes_armatus         -1.3286935          0.8097102
## Parasergestes_vigilax         -1.9939834          1.3460569
## Phorcosergia_grandis          -1.3527442          0.8688288
## Robustosergia_regalis         -1.2356712          0.8429113
## Robustosergia_robusta         -1.1602404          0.8095839
## Sergestes_atlanticus          -1.8546790          1.2481074
## Sergia_tenuiremis             -0.8172119          0.5507995
## 
## attr(,"class")
## [1] "coef.mer"
# AIC
AIC(lme_nopres)
## [1] -1123.278
#fixed effects
fixef(lme_nopres)
##        (Intercept) log10(Body_Length) 
##         -1.3519333          0.8937128
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_nopres, level = 0.95)
##                          2.5 %      97.5 %
## .sig01              0.27582268  0.74974134
## .sig02             -0.99580165 -0.97348804
## .sig03              0.17701128  0.50987568
## .sigma              0.05903897  0.06767603
## (Intercept)        -1.64344107 -1.07983293
## log10(Body_Length)  0.71455335  1.09409513

Model including preservation type

#reorder levels for preservation so that contrasts will be compared to fresh samples
specimens_test$Preservation <- factor(specimens_test$Preservation, 
                                      levels = c("fresh", "ethanol", "paraformaldehyde"))

# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_withpres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + log10(Body_Length) | genus_species), 
                   data = specimens_test,
                   REML = FALSE) #uses ML instead of REML

# summary
summary(lme_withpres)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 +  
##     log10(Body_Length) | genus_species)
##    Data: specimens_test
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1152.9  -1120.0    584.5  -1168.9      442 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7037 -0.5122  0.0295  0.5644  4.3676 
## 
## Random effects:
##  Groups        Name               Variance Std.Dev. Corr 
##  genus_species (Intercept)        0.167876 0.40973       
##                log10(Body_Length) 0.071632 0.26764  -0.99
##  Residual                         0.003701 0.06083       
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                               Estimate Std. Error t value
## (Intercept)                  -1.353719   0.119794 -11.300
## log10(Body_Length)            0.894613   0.078169  11.445
## Preservationethanol          -0.019607   0.008151  -2.405
## Preservationparaformaldehyde  0.045520   0.012167   3.741
## 
## Correlation of Fixed Effects:
##             (Intr) l10(B_ Prsrvtnt
## lg10(Bdy_L) -0.990                
## Prsrvtnthnl -0.057  0.022         
## Prsrvtnprfr -0.068  0.041  0.455
# AIC
AIC(lme_withpres)
## [1] -1152.918
#fixed effects
fixef(lme_withpres)
##                  (Intercept)           log10(Body_Length) 
##                  -1.35371863                   0.89461257 
##          Preservationethanol Preservationparaformaldehyde 
##                  -0.01960702                   0.04551984
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_withpres, level = 0.95)
##                                    2.5 %      97.5 %
## .sig01                        0.25281957  0.67323848
## .sig02                       -0.99341603 -0.96585932
## .sig03                        0.15905144  0.44860338
## .sigma                        0.05692248  0.06521126
## (Intercept)                  -1.61782966 -1.10480110
## log10(Body_Length)            0.73401477  1.07261071
## Preservationethanol          -0.03573800 -0.00350345
## Preservationparaformaldehyde  0.02141721  0.06962976

Compare models

Then we can compare the two models via AIC and log likelihood.

#model AIC comparison
AIC(lme_nopres, lme_withpres)
##              df       AIC
## lme_nopres    6 -1123.278
## lme_withpres  8 -1152.918
#log likelihood comparison
anova(lme_nopres, lme_withpres)
## Data: specimens_test
## Models:
## lme_nopres: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species)
## lme_withpres: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + log10(Body_Length) | genus_species)
##              npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)    
## lme_nopres      6 -1123.3 -1098.6 567.64  -1135.3                        
## lme_withpres    8 -1152.9 -1120.0 584.46  -1168.9 33.64  2  4.957e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here we see that the model that includes preservation type as a covariate is significantly better supported. Ethanol preservation has a negative effect on intercept (eye size compared to body size), but paraformaldehyde has a positive effect. Note also that these effects are very small, so while they are significantly affecting our models, they seem to have minor effects.

Linear SMA regression for pairwise comparisons among species slopes

The smatr package also allows you to test whether groups exhibit significant differences in allometric scaling (slopes or intercepts). Below, we again use standardized major axis regression to model the scaling of eye diameter with body length across all specimens measured and test whether we see significant differences in static allometry across species.

# Model for static allometry across specimens with species and species*body size interaction as covariates
sma_species <- sma(formula = Eye_Diameter ~ Body_Length * genus_species, 
            data = specimens_test, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            multcomp = TRUE, # adds pairwise comparisons between groups
            multcompmethod = "adjusted", # adjusts p-value for multiple comparisons
            alpha = 0.05) 

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_species, which = "default",type = "o") 
plot(sma_species,which = "residual",type = "o")
plot(sma_species, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_species)
## Call: sma(formula = Eye_Diameter ~ Body_Length * genus_species, data = specimens_test, 
##     log = "xy", method = "SMA", alpha = 0.05, multcomp = TRUE, 
##     multcompmethod = "adjusted") 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Results of comparing lines among groups.
## 
## H0 : slopes are equal.
## Likelihood ratio statistic : 78.37 with 14 degrees of freedom
## P-value : 5.6738e-11 
## ------------------------------------------------------------
## 
## Results of multiple comparisons among groups.
## 
## Test for pair-wise difference in slope :
##                  genus_species_1              genus_species_2       Pval
## 1       Allosergestes_pectinatus       Allosergestes_sargassi 1.00000000
## 2       Allosergestes_pectinatus Challengerosergia_hansjacobi 0.10631965
## 3       Allosergestes_pectinatus  Challengerosergia_talismani 0.03771725
## 4       Allosergestes_pectinatus      Deosergestes_corniculum 0.00370375
## 5       Allosergestes_pectinatus         Deosergestes_henseni 0.00191135
## 6       Allosergestes_pectinatus         Eusergestes_arcticus 1.00000000
## 7       Allosergestes_pectinatus    Gardinerosergia_splendens 0.00417512
## 8       Allosergestes_pectinatus        Parasergestes_armatus 0.14662731
## 9       Allosergestes_pectinatus        Parasergestes_vigilax 1.00000000
## 10      Allosergestes_pectinatus         Phorcosergia_grandis 0.00484593
## 11      Allosergestes_pectinatus        Robustosergia_regalis 0.00330606
## 12      Allosergestes_pectinatus        Robustosergia_robusta 0.00386781
## 13      Allosergestes_pectinatus         Sergestes_atlanticus 1.00000000
## 14      Allosergestes_pectinatus            Sergia_tenuiremis 0.00977338
## 15        Allosergestes_sargassi Challengerosergia_hansjacobi 0.17642048
## 16        Allosergestes_sargassi  Challengerosergia_talismani 0.03323729
## 17        Allosergestes_sargassi      Deosergestes_corniculum 0.00219121
## 18        Allosergestes_sargassi         Deosergestes_henseni 0.00062265
## 19        Allosergestes_sargassi         Eusergestes_arcticus 1.00000000
## 20        Allosergestes_sargassi    Gardinerosergia_splendens 0.00279183
## 21        Allosergestes_sargassi        Parasergestes_armatus 0.25381226
## 22        Allosergestes_sargassi        Parasergestes_vigilax 1.00000000
## 23        Allosergestes_sargassi         Phorcosergia_grandis 0.00071873
## 24        Allosergestes_sargassi        Robustosergia_regalis 0.00055799
## 25        Allosergestes_sargassi        Robustosergia_robusta 0.00340588
## 26        Allosergestes_sargassi         Sergestes_atlanticus 1.00000000
## 27        Allosergestes_sargassi            Sergia_tenuiremis 0.01682182
## 28  Challengerosergia_hansjacobi  Challengerosergia_talismani 1.00000000
## 29  Challengerosergia_hansjacobi      Deosergestes_corniculum 1.00000000
## 30  Challengerosergia_hansjacobi         Deosergestes_henseni 0.99999989
## 31  Challengerosergia_hansjacobi         Eusergestes_arcticus 1.00000000
## 32  Challengerosergia_hansjacobi    Gardinerosergia_splendens 1.00000000
## 33  Challengerosergia_hansjacobi        Parasergestes_armatus 1.00000000
## 34  Challengerosergia_hansjacobi        Parasergestes_vigilax 0.04921452
## 35  Challengerosergia_hansjacobi         Phorcosergia_grandis 1.00000000
## 36  Challengerosergia_hansjacobi        Robustosergia_regalis 1.00000000
## 37  Challengerosergia_hansjacobi        Robustosergia_robusta 1.00000000
## 38  Challengerosergia_hansjacobi         Sergestes_atlanticus 0.25203409
## 39  Challengerosergia_hansjacobi            Sergia_tenuiremis 1.00000000
## 40   Challengerosergia_talismani      Deosergestes_corniculum 0.99999999
## 41   Challengerosergia_talismani         Deosergestes_henseni 0.99993757
## 42   Challengerosergia_talismani         Eusergestes_arcticus 1.00000000
## 43   Challengerosergia_talismani    Gardinerosergia_splendens 1.00000000
## 44   Challengerosergia_talismani        Parasergestes_armatus 1.00000000
## 45   Challengerosergia_talismani        Parasergestes_vigilax 0.02124209
## 46   Challengerosergia_talismani         Phorcosergia_grandis 1.00000000
## 47   Challengerosergia_talismani        Robustosergia_regalis 1.00000000
## 48   Challengerosergia_talismani        Robustosergia_robusta 0.99999992
## 49   Challengerosergia_talismani         Sergestes_atlanticus 0.11697344
## 50   Challengerosergia_talismani            Sergia_tenuiremis 0.99999999
## 51       Deosergestes_corniculum         Deosergestes_henseni 1.00000000
## 52       Deosergestes_corniculum         Eusergestes_arcticus 0.99999999
## 53       Deosergestes_corniculum    Gardinerosergia_splendens 1.00000000
## 54       Deosergestes_corniculum        Parasergestes_armatus 1.00000000
## 55       Deosergestes_corniculum        Parasergestes_vigilax 0.00236873
## 56       Deosergestes_corniculum         Phorcosergia_grandis 1.00000000
## 57       Deosergestes_corniculum        Robustosergia_regalis 1.00000000
## 58       Deosergestes_corniculum        Robustosergia_robusta 1.00000000
## 59       Deosergestes_corniculum         Sergestes_atlanticus 0.01350778
## 60       Deosergestes_corniculum            Sergia_tenuiremis 1.00000000
## 61          Deosergestes_henseni         Eusergestes_arcticus 0.99999959
## 62          Deosergestes_henseni    Gardinerosergia_splendens 1.00000000
## 63          Deosergestes_henseni        Parasergestes_armatus 0.99999978
## 64          Deosergestes_henseni        Parasergestes_vigilax 0.00128436
## 65          Deosergestes_henseni         Phorcosergia_grandis 1.00000000
## 66          Deosergestes_henseni        Robustosergia_regalis 1.00000000
## 67          Deosergestes_henseni        Robustosergia_robusta 1.00000000
## 68          Deosergestes_henseni         Sergestes_atlanticus 0.00731827
## 69          Deosergestes_henseni            Sergia_tenuiremis 1.00000000
## 70          Eusergestes_arcticus    Gardinerosergia_splendens 0.99999997
## 71          Eusergestes_arcticus        Parasergestes_armatus 1.00000000
## 72          Eusergestes_arcticus        Parasergestes_vigilax 0.99996770
## 73          Eusergestes_arcticus         Phorcosergia_grandis 1.00000000
## 74          Eusergestes_arcticus        Robustosergia_regalis 1.00000000
## 75          Eusergestes_arcticus        Robustosergia_robusta 0.99999993
## 76          Eusergestes_arcticus         Sergestes_atlanticus 1.00000000
## 77          Eusergestes_arcticus            Sergia_tenuiremis 0.99999911
## 78     Gardinerosergia_splendens        Parasergestes_armatus 1.00000000
## 79     Gardinerosergia_splendens        Parasergestes_vigilax 0.00266011
## 80     Gardinerosergia_splendens         Phorcosergia_grandis 1.00000000
## 81     Gardinerosergia_splendens        Robustosergia_regalis 1.00000000
## 82     Gardinerosergia_splendens        Robustosergia_robusta 1.00000000
## 83     Gardinerosergia_splendens         Sergestes_atlanticus 0.01499340
## 84     Gardinerosergia_splendens            Sergia_tenuiremis 1.00000000
## 85         Parasergestes_armatus        Parasergestes_vigilax 0.06386567
## 86         Parasergestes_armatus         Phorcosergia_grandis 1.00000000
## 87         Parasergestes_armatus        Robustosergia_regalis 1.00000000
## 88         Parasergestes_armatus        Robustosergia_robusta 1.00000000
## 89         Parasergestes_armatus         Sergestes_atlanticus 0.31517941
## 90         Parasergestes_armatus            Sergia_tenuiremis 1.00000000
## 91         Parasergestes_vigilax         Phorcosergia_grandis 0.00320670
## 92         Parasergestes_vigilax        Robustosergia_regalis 0.00220257
## 93         Parasergestes_vigilax        Robustosergia_robusta 0.00240837
## 94         Parasergestes_vigilax         Sergestes_atlanticus 1.00000000
## 95         Parasergestes_vigilax            Sergia_tenuiremis 0.00481337
## 96          Phorcosergia_grandis        Robustosergia_regalis 1.00000000
## 97          Phorcosergia_grandis        Robustosergia_robusta 1.00000000
## 98          Phorcosergia_grandis         Sergestes_atlanticus 0.01851043
## 99          Phorcosergia_grandis            Sergia_tenuiremis 1.00000000
## 100        Robustosergia_regalis        Robustosergia_robusta 1.00000000
## 101        Robustosergia_regalis         Sergestes_atlanticus 0.01271145
## 102        Robustosergia_regalis            Sergia_tenuiremis 1.00000000
## 103        Robustosergia_robusta         Sergestes_atlanticus 0.01358845
## 104        Robustosergia_robusta            Sergia_tenuiremis 1.00000000
## 105         Sergestes_atlanticus            Sergia_tenuiremis 0.02662028
##         TestStat
## 1   7.666232e-02
## 2   1.070239e+01
## 3   1.269776e+01
## 4   1.710663e+01
## 5   1.836679e+01
## 6   1.514543e+00
## 7   1.687873e+01
## 8   1.006766e+01
## 9   8.341304e-01
## 10  1.659541e+01
## 11  1.732280e+01
## 12  1.702417e+01
## 13  8.875795e-04
## 14  1.526295e+01
## 15  9.695878e+00
## 16  1.293861e+01
## 17  1.810626e+01
## 18  2.051006e+01
## 19  1.218474e+00
## 20  1.764469e+01
## 21  8.943571e+00
## 22  1.268181e+00
## 23  2.023537e+01
## 24  2.072006e+01
## 25  1.726619e+01
## 26  7.817618e-02
## 27  1.423248e+01
## 28  2.850713e-03
## 29  1.459201e+00
## 30  2.159382e+00
## 31  7.165782e-01
## 32  1.474497e+00
## 33  3.154991e-03
## 34  1.218978e+01
## 35  7.281210e-01
## 36  1.215204e+00
## 37  1.670749e+00
## 38  8.958462e+00
## 39  1.652888e+00
## 40  1.971056e+00
## 41  2.909088e+00
## 42  7.908634e-01
## 43  1.901819e+00
## 44  1.381453e-02
## 45  1.378961e+01
## 46  1.203054e+00
## 47  1.877144e+00
## 48  2.130775e+00
## 49  1.051505e+01
## 50  1.917202e+00
## 51  1.091006e-01
## 52  1.922892e+00
## 53  9.383239e-03
## 54  1.528670e+00
## 55  1.795778e+01
## 56  4.919054e-01
## 57  9.043887e-02
## 58  3.949879e-02
## 59  1.464885e+01
## 60  1.486450e-01
## 61  2.285326e+00
## 62  4.454856e-02
## 63  2.224689e+00
## 64  1.912552e+01
## 65  1.138002e+00
## 66  4.397057e-01
## 67  1.181780e-02
## 68  1.581216e+01
## 69  1.724888e-02
## 70  2.036867e+00
## 71  6.530716e-01
## 72  2.808168e+00
## 73  1.445758e+00
## 74  1.716984e+00
## 75  2.123385e+00
## 76  1.501074e+00
## 77  2.362523e+00
## 78  1.546238e+00
## 79  1.773675e+01
## 80  5.367445e-01
## 81  1.427310e-01
## 82  9.239058e-03
## 83  1.445085e+01
## 84  8.641470e-02
## 85  1.169009e+01
## 86  7.995621e-01
## 87  1.285411e+00
## 88  1.740929e+00
## 89  8.475706e+00
## 90  1.723105e+00
## 91  1.738089e+01
## 92  1.809640e+01
## 93  1.792616e+01
## 94  6.837798e-01
## 95  1.660823e+01
## 96  2.447817e-01
## 97  7.055402e-01
## 98  1.405093e+01
## 99  7.801888e-01
## 100 2.376709e-01
## 101 1.476415e+01
## 102 3.827565e-01
## 103 1.463755e+01
## 104 4.529580e-02
## 105 1.336087e+01
## 
## ------------------------------------------------------------
## Coefficients by group in variable "genus_species"
## 
## Group: Allosergestes_pectinatus 
##             elevation    slope
## estimate    -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
## 
## H0 : variables uncorrelated.
## R-squared : 0.8156794 
## P-value : 9.643e-06 
## 
## Group: Allosergestes_sargassi 
##             elevation    slope
## estimate    -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
## 
## H0 : variables uncorrelated.
## R-squared : 0.4037707 
## P-value : 2.3852e-05 
## 
## Group: Challengerosergia_hansjacobi 
##             elevation    slope
## estimate    -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
## 
## H0 : variables uncorrelated.
## R-squared : 0.4129122 
## P-value : 3.3885e-06 
## 
## Group: Challengerosergia_talismani 
##             elevation     slope
## estimate    -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
## 
## H0 : variables uncorrelated.
## R-squared : 0.8120894 
## P-value : 2.0471e-12 
## 
## Group: Deosergestes_corniculum 
##             elevation     slope
## estimate    -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
## 
## H0 : variables uncorrelated.
## R-squared : 0.8421211 
## P-value : 8.1768e-08 
## 
## Group: Deosergestes_henseni 
##             elevation     slope
## estimate    -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
## 
## H0 : variables uncorrelated.
## R-squared : 0.3380975 
## P-value : 1.6982e-06 
## 
## Group: Eusergestes_arcticus 
##              elevation     slope
## estimate    -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
## 
## H0 : variables uncorrelated.
## R-squared : 0.8368056 
## P-value : 0.029484 
## 
## Group: Gardinerosergia_splendens 
##              elevation     slope
## estimate    -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
## 
## H0 : variables uncorrelated.
## R-squared : 0.4597308 
## P-value : 1.1132e-06 
## 
## Group: Parasergestes_armatus 
##             elevation     slope
## estimate    -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
## 
## H0 : variables uncorrelated.
## R-squared : 0.5670657 
## P-value : 1.5744e-06 
## 
## Group: Parasergestes_vigilax 
##             elevation    slope
## estimate    -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
## 
## H0 : variables uncorrelated.
## R-squared : 0.5274941 
## P-value : 0.00096145 
## 
## Group: Phorcosergia_grandis 
##             elevation     slope
## estimate    -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
## 
## H0 : variables uncorrelated.
## R-squared : 0.8798127 
## P-value : < 2.22e-16 
## 
## Group: Robustosergia_regalis 
##             elevation     slope
## estimate    -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
## 
## H0 : variables uncorrelated.
## R-squared : 0.8280327 
## P-value : 2.4971e-15 
## 
## Group: Robustosergia_robusta 
##              elevation     slope
## estimate    -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
## 
## H0 : variables uncorrelated.
## R-squared : 0.6989428 
## P-value : 6.7833e-07 
## 
## Group: Sergestes_atlanticus 
##             elevation    slope
## estimate    -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
## 
## H0 : variables uncorrelated.
## R-squared : 0.7067609 
## P-value : 8.6609e-05 
## 
## Group: Sergia_tenuiremis 
##              elevation     slope
## estimate    -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
## 
## H0 : variables uncorrelated.
## R-squared : 0.2880646 
## P-value : 0.0032339
#save coefficients of fits as object
cc_species <- data.frame(coef(sma_species))

Among 15 sergestid species (450 specimens), species significantly differ in eye-body allometric scaling (p < 0.0001). Pairwise comparisons show which species differ significantly in slope and which do not.

Here we plot the static allometry of all the species with sufficient sampling (n = 15). It’s interactive, so you can click on different species in the legend to hide or show those points.

#reorder factor levels for figure legend (phylo order)
specimens_test$species_text <- factor(specimens_test$species_text, 
                                      levels = c("Deosergestes corniculum",
                                                 "Deosergestes henseni",
                                                 "Allosergestes pectinatus",
                                                 "Allosergestes sargassi",
                                                 "Sergestes atlanticus",
                                                 "Parasergestes vigilax",
                                                 "Parasergestes armatus",
                                                 "Eusergestes arcticus",
                                                 "Gardinerosergia splendens",
                                                 "Robustosergia regalis",
                                                 "Robustosergia robusta",
                                                 "Phorcosergia grandis",
                                                 "Sergia tenuiremis",
                                                 "Challengerosergia talismani",
                                                 "Challengerosergia hansjacobi"))

#make shape pallette
shapes.sp <- c("Deosergestes corniculum" = 21,
              "Deosergestes henseni" = 22, 
              "Allosergestes pectinatus" = 23, 
              "Allosergestes sargassi" = 24,
              "Sergestes atlanticus" = 25, 
              "Parasergestes vigilax" = 21, 
              "Parasergestes armatus" = 22,
              "Eusergestes arcticus" = 23, 
              "Gardinerosergia splendens" = 24, 
              "Robustosergia regalis" = 25, 
              "Robustosergia robusta" = 21, 
              "Phorcosergia grandis" = 22,
              "Sergia tenuiremis" = 23,
              "Challengerosergia talismani" = 24,
              "Challengerosergia hansjacobi" = 25)

#rainbow palette
cols.sp <- c("Deosergestes corniculum" = "#16ABA8",
              "Deosergestes henseni" = "#1f78b4",
              "Allosergestes pectinatus" = "#3acf75",
              "Allosergestes sargassi" = "#33a02c",
              "Sergestes atlanticus" = "#fb9a99",
              "Parasergestes vigilax" = "#e31a1c",
              "Parasergestes armatus" = "#fdbf6f",
              "Eusergestes arcticus" = "#ff7f00",
              "Gardinerosergia splendens" = "#cab2d6",
              "Robustosergia regalis" = "#6a3d9a",
              "Robustosergia robusta" = "gray31",
              "Phorcosergia grandis" = "burlywood4",
              "Sergia tenuiremis" = "black",
              "Challengerosergia talismani" = "maroon1",
              "Challengerosergia hansjacobi" = "orchid1")

#sergia/sergestes green purple pallette
# cols.sp <- c(#Sergestes group
#               "Deosergestes corniculum" = "#512E5F",
#               "Deosergestes henseni" = "#633974", 
#               "Allosergestes pectinatus" = "#76448A", 
#               "Allosergestes sargassi" = "#884EA0",
#               "Sergestes atlanticus" = "#9B59B6", 
#               #"Neosergestes edwardsii" = "#AF7AC5"
#               "Parasergestes vigilax" = "#C39BD3", 
#               "Parasergestes armatus" = "#D7BDE2",
#               "Eusergestes arcticus" = "#EBDEF0", 
#              #Sergia group
#               "Gardinerosergia splendens" = "#D4EFDF",
#               "Robustosergia regalis" = "#A9DFBF", 
#               "Robustosergia robusta" = "#52BE80",
#               "Phorcosergia grandis" = "#27AE60",
#               "Sergia tenuiremis" = "#1E8449",
#               "Challengerosergia talismani" = "#196F3D",
#               "Challengerosergia hansjacobi" = "#145A32")


#set s limits by species for plots
limits <- specimens_test %>%
  group_by(genus_species) %>%
  summarise(max = max(Body_Length, na.rm = TRUE),
            min = min(Body_Length, na.rm = TRUE))

#make plot
plot_species <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color =  species_text, shape = species_text, fill = species_text)) + 
  geom_point(size = 2, alpha = 0.4) + 
  scale_shape_manual(values = shapes.sp, name = "Species") + 
  scale_color_manual(values = cols.sp, name = "Species") +
  scale_fill_manual(values = cols.sp, name = "Species") +
  scale_x_log10(name = "Body length (mm)", breaks = c(15,20,30,50,100)) + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)", breaks = c(0.2,0.3,0.5,1, 2.6)) + #makes y axis log scale and named
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.text = element_text(face = "italic")) +
    geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_pectinatus")], xend = limits$max[which(limits$genus_species == "Allosergestes_pectinatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"])), color = cols.sp["Allosergestes pectinatus"]) + #Rana temporaria adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_sargassi")], xend = limits$max[which(limits$genus_species == "Allosergestes_sargassi")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"])), color = cols.sp["Allosergestes sargassi"]) + #Allosergestes_sargassi adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")], xend = limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"])), color = cols.sp["Challengerosergia hansjacobi"]) + #Challengerosergia_hansjacobi adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_talismani")], xend = limits$max[which(limits$genus_species == "Challengerosergia_talismani")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"])), color = cols.sp["Challengerosergia talismani"]) + #Challengerosergia_talismani adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_corniculum")], xend = limits$max[which(limits$genus_species == "Deosergestes_corniculum")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"])), color = cols.sp["Deosergestes corniculum"]) + #Deosergestes_corniculum adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_henseni")], xend = limits$max[which(limits$genus_species == "Deosergestes_henseni")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"])), color = cols.sp["Deosergestes henseni"]) + #Deosergestes henseni adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Gardinerosergia_splendens")], xend = limits$max[which(limits$genus_species == "Gardinerosergia_splendens")], y = 10^(log10(limits$min[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"])), color = cols.sp["Gardinerosergia splendens"]) + #Gardinerosergia splendens adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_armatus")], xend = limits$max[which(limits$genus_species == "Parasergestes_armatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"])), color = cols.sp["Parasergestes armatus"]) + #Parasergestes armatus adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_vigilax")], xend = limits$max[which(limits$genus_species == "Parasergestes_vigilax")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"])), color = cols.sp["Parasergestes vigilax"]) + #Parasergestes vigilax adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Phorcosergia_grandis")], xend = limits$max[which(limits$genus_species == "Phorcosergia_grandis")], y = 10^(log10(limits$min[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"])), color = cols.sp["Phorcosergia grandis"]) + #Phorcosergia grandis adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_robusta")], xend = limits$max[which(limits$genus_species == "Robustosergia_robusta")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_robusta")])*cc_species["Robustosergia_robusta", "slope"] + cc_species["Robustosergia_robusta", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_robusta")])*cc_species["Robustosergia_robusta", "slope"] + cc_species["Robustosergia_robusta", "elevation"])), color = cols.sp["Robustosergia robusta"]) + #Robustosergia robusta adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_regalis")], xend = limits$max[which(limits$genus_species == "Robustosergia_regalis")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"])), color = cols.sp["Robustosergia regalis"]) +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergestes_atlanticus")], xend = limits$max[which(limits$genus_species == "Sergestes_atlanticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"])), color = cols.sp["Sergestes atlanticus"]) +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergia_tenuiremis")], xend = limits$max[which(limits$genus_species == "Sergia_tenuiremis")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"])), color = cols.sp["Sergia tenuiremis"])+
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Eusergestes_arcticus")], xend = limits$max[which(limits$genus_species == "Eusergestes_arcticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"])), color = cols.sp["Eusergestes arcticus"])+
  theme(legend.text = element_text(face = "italic"))


#interactive plot
ggplotly(plot_species)
#make plot with alpha = 1 for legend
plot_species_leg <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color =  species_text, shape = species_text, fill = species_text)) + 
  geom_point(size = 2, alpha = 1) + 
  scale_shape_manual(values = shapes.sp, name = "Species") + 
  scale_color_manual(values = cols.sp, name = "Species") +
  scale_fill_manual(values = cols.sp, name = "Species") +
  theme_bw() +
  theme(legend.text = element_text(face = "italic"))

# remove legend from previous plot
plot <- plot_species + theme(legend.position = "none")

# extract legend from this plot with alpha = 1
leg <- get_legend(plot_species_leg)

#export figure
#png("../Figures/dev-species.png", width = 12, height = 12, units = "in", res = 500)
pdf("../Figures/Fig-1.pdf", width = 9, height = 5)
plot_grid(plot, leg, ncol = 2, rel_widths = c(1, .4))
dev.off()